
#ifndef VECTOR_HPP
#define VECTOR_HPP

// $Id: Vector.hpp 2 2012-05-22 17:55:16Z gerry@gaboury.biz $

#include <vector>
#include <algorithm>
#include <stdexcept>
#include <iostream>
#include <iomanip>
//#include "VectorStore.hpp"

using namespace std;

template <class T> class Vector : public vector<T>
{
  using vector<T>::iterator;
  //using vector<T>::resize;

public:

  Vector () : vector<T> () { }

  Vector ( int p_size ) : vector<T> ( p_size ) { }

  // initialize with a value
  Vector ( int p_size, T p_val ) : vector<T> ( p_size, p_val ) { }

  // initialize with an array
  Vector ( int p_size, T* p_arr ) : vector<T> ( p_size ) 
  {
    for (int i=0; i < p_size; ++i )
      (*this)[i] = p_arr[i];
  }

  Vector<T> operator+ ( const Vector<T>& p_vec )
  {
    if ( p_vec.size() < this->size() )
      throw (range_error("Vector::operator+ : range error"));

    // the storage class
    //VectorStore* VS = VectorStore::instance();
    //Vector<T>* Vtmp = VS->newVector<T>( this->size() );

    Vector<T> Vtmp = Vector<T>( this->size() );

    for (unsigned j = 0; j < this->size(); ++j )
      {
	Vtmp[j] = (*this)[j] + p_vec[j];
      }
    return Vtmp;
  }

  Vector<T> operator-(const Vector<T>& p_vec)
  {
    if ( p_vec.size() < this->size() )
      throw (range_error("Vector::operator- : range error"));

    // the storage class
    //VectorStore* VS = VectorStore::instance();
    //Vector<T>* Vtmp = VS->newVector<T>( this->size() );
    Vector<T> Vtmp ( this->size() );

    for (unsigned j = 0; j < this->size(); ++j )
      {
	Vtmp[j] = (*this)[j] - p_vec[j];
      }
    return Vtmp;
  }

  T operator*(const Vector<T>& p_vec)
  {
    if ( p_vec.size() < this->size() )
      throw (range_error("Vector::operator* : range error"));

    T tmp = 0;
    for ( unsigned j = 0; j < this->size(); ++j )
      {
	tmp += (*this)[j] * p_vec[j];
      }
    return tmp;
  }

  Vector<T>& operator= ( const Vector& p_vec )
  {
    if ( this->size() != p_vec.size() )
      resize ( p_vec.size() );

    for (unsigned j = 0; j < p_vec.size(); ++j )
      {
	(*this)[j] = p_vec[j];
      }
    return *this;
  }

  Vector<T>& operator+= ( const Vector& p_vec )
  {
    if ( this->size() != p_vec.size() )
      throw (range_error("Vector::operator+= : size mismatch"));

    for (unsigned j = 0; j < p_vec.size(); ++j )
      {
	(*this)[j] += p_vec[j];
      }
    return *this;
  }

  Vector<T>& operator-= ( const Vector& p_vec )
  {
    if ( this->size() != p_vec.size() )
      throw (range_error("Vector::operator+= : size mismatch"));

    for (unsigned j = 0; j < p_vec.size(); ++j )
      {
	(*this)[j] -= p_vec[j];
      }
    return *this;
  }

  Vector<T> operator* ( const T& p_scalar )
  {
    Vector<T> Vtmp ( this->size() );
    for (unsigned j = 0; j < this->size(); ++j )
      {
	Vtmp[j] = (*this)[j] * p_scalar;
      }
    return Vtmp;
  }

  Vector<T>& operator*= ( const T& p_scalar )
  {
    for (unsigned j = 0; j < this->size(); ++j )
      {
	(*this)[j] *= p_scalar;
      }
    return *this;
  }
  Vector<T>& operator/= ( const T& p_scalar )
  {
    for (unsigned j = 0; j < this->size(); ++j )
      {
	(*this)[j] /= p_scalar;
      }
    return *this;
  }

  Vector<T>& operator= ( const T& p_scalar )
  {
    fill ( this->begin(), this->end(), p_scalar );
    return *this;
  }

  T SumSq() const
  {
    T sum = 0;
    typename vector<T>::const_iterator iter;
    for ( iter = this->begin(); iter != this->end(); ++iter )
      sum += (*iter) * (*iter);
    return sum;
  }
  //
  // mean and sample variance
  //
  pair<T,T> MeanVariance() const
  {
    T mean=0, var=0;
    unsigned n = this->size();
    typename vector<T>::const_iterator iter;
    for ( iter = this->begin(); iter != this->end(); ++iter )
      {
	mean += *iter;
	var += (*iter) * (*iter);
      }
    mean /= (T)n;
    var -= mean * mean * n;
    var /= (T)(n-1);
    return pair<T,T>(mean,var);
  }

  // IO
  void write(ostream *os, int p_prec=5, int p_width=8) const
  {
    // set output format
    ios::fmtflags oldopt = os->setf( ios::fixed );
    os->precision(p_prec);

    for (int i=0; i < this->size(); ++i)
      *os << setw(p_width) << (*this)[i] << " ";
    *os << endl;

    // reset output format
    os->flags( oldopt );
  }
  void read(istream *is)
  {
    for (int i=0; i < this->size(); ++i)
      *is >> (*this)[i];
  }
private:
  // disable these methods
  //virtual void push_back (const T&) { }
  virtual void pop_back () { }
};

template <class T> ostream &operator<< ( ostream &os,Vector<T> &vec )
{
  vec.write(&os);
  return os;
}
template <class T> ostream &operator<< ( ostream &os,const Vector<T> &vec )
{
  vec.write(&os);
  return os;
}
template <class T> istream &operator>> ( istream &is, Vector<T> &vec )
{
  vec.read(&is);
  return is;
}

#endif // VECTOR_HPP
